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We present a method for bridging the gap between the Dirac effective field theory and atom- 
istic simulations in graphene based on the Husimi projection, allowing us to depict phenomena in 
graphene at arbitrary scales. This technique takes the atomistic wavef unction as an input, and 
produces semiclassical pictures of quasiparticles in the two Dirac valleys. We use the Husimi tech- 
nique to produce maps of the scattering behavior of boundaries, giving insight into the properties 
of wavefunctions at energies both close to and far from the Dirac point. Boundary conditions play 
a significant role to the rise of Fano resonances, which we examine using the Husimi map to deepen 
our understanding of bond currents near resonance. 



I. INTRODUCTION 

With interest and experimental capabilities in 
graphene devices growing[l -8j, the need has never been 
greater to improve our understanding of quantum states 
in this material. Despite the success of the Dirac effective 
field theory for graphene[9j, however, many technological 
proposals arise from predictions using the more funda- 
mental tight-binding approximation [10- 13J. This is be- 
cause the atomistic model that underlies the Dirac theory 
is able to incorporate phenomena such as scattering from 
small defects [I3HI51 . ripplespjj], or edge types [2QH22] - 
all of which promise technological applications. How- 
ever, atomistic calculations are computationally expen- 
sive, and replacing these features with scattering theo- 
ries in a more-efficient Dirac model introduces substan- 
tial challenges. A robust approach which can analyze the 
atomistic wavefunction to produce semiclassical pictures 
of quasiparticles in the two Dirac valleys remains to be 
seen. 

To address these issues and expand our understanding 
of graphene quantum states, we use the Husimi projec- 
tion technique, introduced by Mason et al. [23— 25J , to 
produce snapshots of the local momentum distribution 
and underlying semiclassical structure in graphene wave- 
functions. When Husimi projections are calculated at 
many points across a system, the Husimi map that results 
provides a semiclassical picture of the atomistic wave- 
function. In this article, we define the Husimi map for 
graphene systems (Sec.|n|), and use it to deepen our un- 
derstanding of boundary conditio ns in b oth high-energy 
relativistic scar state s[26, 2 7| (Sec. Ill A), and states near 
the Dirac point (Sec. IIIB). We then use Husimi maps to 
interpret Fano resonances[28-30| within this novel mate- 
rial (Sec. gncj. 



II. METHOD 

A. Definition of the Husimi Projection 

The conduction band of the graphene system can be 
approximated as a honeycomb lattice with a single p z 
orbital located at each carbon- atom lattice site 1 9 . The 



Husimi function is defined as a measurement between a 
wavefunction ^({r^}) defined at each orbital, and a co- 
herent state |ro,ko,cr) which describes an envelope func- 
tion over those sites that minimizes the joint uncertainty 
in spatial and momentum coordinates. The parameter <j 
defines the spatial spread of the coherent state and de- 
fines the uncertainties in space and momentum according 
to the well-known relation 
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As a result, there is a trade-off for any value of a selected: 
for small a, there is better spatial resolution but poorer 
resolution in /c-space, and vice versa for large a. 

Writing out the dot product of the wavefunction and 
the coherent state 



(^|r ,k ,cr) 



the Husimi function is defined as 
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Hu(r ,k ,cr;^({ri})) = |(^| r , k , cr)\ 2 . (3) 

Weighting the Husimi function by the wavevector ko 
produces the /c-space Husimi vector, and weighting it by 
the group velocity vector Vk-E (k') produces the group- 
velocity Husimi vector. The latter is a stronger reflection 
of classical dynamics in the system, and is used for all 
results in this paper. At each point in the system, we 
can sweep through /c-space by rotating the wavevector ko 
along the Fermi surface in the dispersion relation. The 
multiple Husimi vectors which result form the full Husimi 
projection, providing a snapshot of the local momentum 
distribution. This paper uses 32 wavevectors along the 
Fermi surface of two-dimensional graphene to produce 
group- velocity Husimi project ions [25], 

Even though a few plane waves may dominate the 
wavefunction, momentum uncertainty of the coherent 
state can result in many non- vanishing Husimi vectors. 
Assuming that the dominant plane waves at a point are 
sufficiently separated in /c-space, it is possible to recover 
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their wavevectors using the Multi-Modal Algorithm in 
Mason et al[25\. This method singles out the impor- 
tant wavevectors contributing to a wavefunction at each 
point; in this paper, we additionally remove results below 
a certain threshold to clarify our results. 

The integral over Husimi vectors at a single point 
defines a new vector-valued function Hu (i-q, <r; ^({r})), 
which is equal to 



Hu (r 



>cr;^({ r J)) = J |(^|r ,ko,cr)| 2 k 



d d k . (4) 



It has been shown that for dfc « 1, this function is equal 
to the flux operator [24J. To better represent the classi- 
cal dynamics of the system we can instead weight the 
integrand by the group velocity Vk^ (k') to obtain the 
group- velocity Husimi flux Hu g (i*o, <r; ^(r)) equal to 



Hu g (r ,a; V>({rJ)) 
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r ,ko,a)| 2 V k £(k )^o. 
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which is used throughout the paper. 

Even though the Husimi projection is related to the 
flux operator, it provides much more information since it 
can be used on stationary states which exhibit zero flux, 
and because it can isolate individual bands and valleys in 
the dispersion relation. Sampling the Husimi projection 
at many points across a system to produce a Husimi map, 
we can produce a much better picture of the classical 
dynamics underlying the wavefunction. 
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Figure 1: A magnified view of a boundary on a graphene 
flake. The orientation of the cut relative to the orientation of 
the lattice can produce two edge types, zigzag (highlighted in 
blue) and armchair (highlighted in red). The two sublattices 
of the unit cell are indicated in black (A-sublattice) and grey 
(B-sublattice). 
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B. The Honeycomb Band Structure 

This paper examines the honeycomb lattice Hamilto- 
nian using the nearest-neighbor tight-binding approxima- 
tion 
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where a| is the creation operator at orbital site z, and 
we sum over the set of nearest neighbors. To compare 
against experiment, the hopping integral value is given by 
t = 2.7eV, while e is set to the value of the Fermi energy [1, 
9J. Eigenstates of closed stadium billiard systems are 
computed using sparse matrix eigensolvers to produce 
individual wavefunctions. 

We study finite graphene systems extracted from an 
infinite honeycomb lattice. A filter is applied to remove 
atom sites which are attached to only one other atom site, 
and to bridge under-coordinated sites whose tt orbitals 
would strongly overlap. As a result, each edge is either a 
pure zig-zag, armchair, or mixed boundary, as shown in 
Fig. [T] Recent studies have suggested that under certain 
circumstances, zig-zag edges reconstruct to form a 5-7 
chain[31j, however their scattering properties appear to 
be identical to regular zig-zag boundaries [32]. We have 
elected not to incorporate these features and leave them 
to future work. 



Figure 2: The two-dimensional dispersion relation for 
graphene demonstrates the two inequivalent valleys as cones 
where the edges of the Brillouin Zones (black lines) meet. 
Dashed white lines indicate the one-dimensional dispersion 
surface at E — 0.5t, while solid white lines indicate E — 0.98£, 
demonstrating extreme triagonal warping. 



The band structure for graphene prominently features 
the two inequivalent K' and K valleys in the energy range 
of — t < E < t[9j, as can be seen in Fig. [2] At energies 
close to the Dirac point E = 0, these valleys exhibit a 
linear dispersion relation and the electron behaves as a 
four-component spinor Dirac particle (two pseudo-spins, 
and two traditional spins). Using the creation operators 
and on the A- and B-sublattices respectively (see 
Fig. [I]) , the two pseudospinors can be written as 

V>±, K (k) = -^=(e-^/ 2 at±e^/ 2 bt) (7) 
tf ±iK <(k) = -^(e^/Vie-^V), (8) 

where 0^ = arctan ^| 2L ^,q = k — and the db signs in- 
dicate whether the positive- or negative-energy solutions 
are being used[9j. While the linear dispersion no longer 
applies at energies above ~ 0.4t, the Dirac basis remains 
useful as a means of describing the classical dynamics 
of graphene throughout the energy range —t<E<t. 
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States near the Dirac point and at the upper edge of this 
spectrum are examined in this paper. 

It might be tempting to obtain a representation of ei- 
ther valley in a graphene wavefunction by subtracting off 
a plane wave whose wavevector corresponds to the ori- 
gin of either K or K' valley, leaving behind the residual 
q = k — K^. However, this approach only works when 
quasiparticles are present in only one valley, an assump- 
tion that cannot be generally guaranteed. 

On the other hand, since wavevectors for each valley 
are sufficiently separated in /c-space, the Husimi projec- 
tion can distinguish each valley unambiguously for most 
momentum uncertainties. Because the valleys are part 
of the same band, a scattered quasiparticle from one val- 
ley can emerge in the other [33J. When this occurs, the 
Husimi map shows quasiparticles in one valley funneling 
into a drain, and quasiparticles in the other valley emit- 
ting from a source at the same point, leaving behind a 
signature for inter-valley scattering. 

Between —t<E<t, the Fermi energy contours warp 
from a circular shape near the Dirac point to trigonal 
contours, which emphasize three directions for each valley 
in the distribution of group velocities v g = A^E (k). As 
a result, the magnitude of the wavevector q = k — 
depends on its orientation: It is bounded above by 
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and from below by 
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(9) 
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When characterizing the momentum uncertainty, we use 
the average of these two quantities. 



III. RESULTS 




Figure 3: The Husimi map for three eigenstates of the closed 
graphene stadium billiard with 20270 orbital sites at energy 
E = 0.974t(a),0.964t(b), and 0.951t(c). All three calcula- 
tions use coherent states with relative uncertainty Ak/k = 
30%, whose breadth is indicated by the double arrows on the 
right. Only the upper-right quarter of each stadia is shown. 
At left, the multi-modal analysis for the K' valley, and at 
right the wavefunction representation. The divergence of the 
Husimi map is indicated in green (red) to for positive (nega- 
tive) values. Angular deflection is indicated in blue (Eq. 12). 
Red boxes indicate the magnified views in Fig. [4] 



A. States Away from the Dirac Point 

Fig. [3] shows Husimi maps for three eigen- 
states of a large closed-system stadium bil- 
liard with 20270 orbital sites at energies of 
E = 0.974t(a),0.964t(b), and 0.951t(c). We have 
chosen these states because they exhibit very clear linear 
trajectories. At energies close to E = t, trajectories ex- 
hibit pronounced trigonal warping, as seen by the three 
preferred directions. While the classical trajectories 
are obvious in the wavefunction itself, the Husimi map 
identifies the direction of each trajectory with respect to 
each valley. 

The presence a few dominant classical paths in each 
wavefunction in Fig. [3] allows us to infer the relationship 
between boundary types and scattering among the two 
Dirac valleys. When a quasiparticle in one valley scatters 
into the other, it appears in the Husimi map as drain. 



We can measure this by summing the divergence for all 
angles in the Husimi map as 

Q div . (r; *) = J D(r,k;<Z>)\V k E(k)\d d k', (11) 

where Z}(r,k;\I/) is defined as the divergence of the 
Husimi map for one wavevector k, 



£>(r,k; *) 



£ 

i=i 



Hu (k, r'^j-Hu (k,r;tf) 



(r' - r) • e* 



x exp 



(r'-rf 
2a 2 
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(12) 



where we sum over the d orthogonal dimensions each as- 
sociated with unit vector e^. The divergence in the K' 
valley, seen in green and red (for positive and negative 
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Figure 4: Magnified views of the divergence and angular de- 
flection in Fig. [3] (red boxes). The sources and drains in 
the X'-valley Husimi map are actually inter-valley scatter- 
ing points, which occur along non zig-zag boundaries. In 
contrast, points of angular deflection that are not sources or 
drains correspond to intra-\a\\ey scatterers and occur along 
pure or nearly-pure zig-zag boundaries. 



values, respectively) in Figs. [3] and [4j shows that the scat- 
tering points all lie along non-zig-zag boundaries. Plots 
for the K valley (not shown) are inverted, corroborating 
the time-reversal symmetry relationship between the two 
valleys. 

On the other hand, when a quasiparticle in one valley 
reflects off a boundary but does not scatter into the other 
valley, the divergence is zero, but the reflection can still 
be measured in the angular deflection of the Husimi map, 



Q ang . (r;*) = j |D abs .(r,k;*)V k £(k)| d d k. 



(13) 



D a bs.( r ? k; S&) is defined as the absolute divergence of the 
Husimi function for one particular trajectory angle with 
a wave vector k, 
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As a result, boundary points with large angular deflection 
are either inter- valley or intra- valley scatterers depending 
on the magnitude of divergence at each point. 

In Figs. [3] and [4j we plot the angular deflection in blue 
to compare to the divergence in green and red. Using this 
information, we can determine that for the wavefunction 
in Fig. |3^l, all boundary scattering points are inter- valley 
scatterers, since all points of angular deflection exhibit 
divergence. The wavefunction in Fig. [3)3, on the other 
hand, only exhibits divergence along the vertical sides 
of the stadium billiard: the horizontal top edge exhibits 
strong angular deflection but no divergence, and consti- 
tutes an intra- valley scatterer. Examining the magnified 
views in Figs. |4^i and [4]3, we see that inter- valley scat- 
terers correspond to armchair edges, and the intra- valley 




Figure 5: A closed-system eigenstate at E — 0.72t for the 
smaller graphene stadium. At top, the filtered Multi-Modal 
analysis is with relative momentum uncertainty Ak/k — 30% 
along with the wavefunction (right). The spread of the co- 
herent state is indicated by the double arrows. At bottom, 
higher-resolution calculations of the divergence (green for pos- 
itive, red for negative) and the angular deflection (blue) are 
shown against the graphene structure. The black circle indi- 
cates where the system boundary is perturbed in the original 



paper[27 as discussed in Section IIIC 



scatterers belong to zig-zag edges, corroborating the find- 
ings at the Dirac point by Akhmerov and Beenakker[34j. 
Similar points of scattering can also be found in Figs. [3J3 
and St. 

Because of the time-reversal relationship between the 
two valleys, the severe restriction on group velocities, and 
the placement of zig-zag and armchair boundaries, no 
path at these energies exists without interacting with an 
inter- valley scatterer (data not shown). By comparison, 
it is not only possible but common to find states near the 
Dirac point that exhibit the opposite: all boundary con- 
ditions which are e xpres sed belong to only m£ra-valley 
scatterers (See Sec. IIIB| ). 

In comparison to Fig. |3j the eigenstate of the much 
smaller graphene stadium system in Fig. [5] does not ap- 
pear to show isolated trajectories in its wavefunction rep- 
resentation. This is not surprising since this system can 
only accommodate five deBroglie wavelengths vertically, 
and three horizontally, severely restricting its ability to 
resolve such trajectories. However, clear self-retracing 
trajectories are quite visible in the Husimi map in Figs.[5j 
with evident sources and drains inhabiting the boundary, 
showing that the Husimi map can yield a semiclassical 
interpretation of the dynamics of the states not possible 
from just the wavefunction of the system. Moreover, be- 
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Figure 6: Schematic indicating the locations of armchair 
(blue) and zig-zag (red) edges in the circular system (left) 
and the Wimmer system (right). 



cause the paths indicated by the Husimi map marshal 
the electron away from lateral boundaries, where leads 
connect to produce the open system in Sec. |III C[ the 
Husimi map helps us understand the role this state plays 
in forming a long-lived resonance in the open system. 

In both Figs. [3] and [5j wavefunctions in graphene away 
from the Dirac point are linked to valley switching classi- 
cal ray paths which bounce back and forth along straight 
lines. These wavefunction enhancements are not strictly 
scars [26J, as first suggested by Huang et a/.[27j, since 
scars are generated by unstable classical periodic orbits in 
the analogous classical limit (group velocity) system. In- 
stead, the wavefunction structures are more likely normal 
quantum confinement to stable zones in classical phase 
space constrained by group- velocity warping at these en- 
ergies. 



B. States Near the Dirac Point 

We now explore the properties of low-energy closed- 
system states in graphene, using the circular graphene 
flake and the distorted circular flake introduced by Wim- 
mer et al.[35\. The latter geometry was chosen because 
its dynamics are chaotic and sensitive to the placement of 
armchair and zig-zag boundaries, which shift as a result 
of the distortion. We indicate the two boundary types 
for both geometries in Fig. [6j. 

In the continuous system, the Fermi wave vector grows 
with the square-root of the energy, but in graphene, the 
effective wavevector q = k — grows linearly. As 
a result, the deBroglie wavelength is much larger for 
the graphene system than for the continuous system at 
similar energy scales, making it difficult to conduct cal- 
culations with sufficient structure in the wavefunction. 
Consequently, we examine states at energies away from 
the Dirac point to bring calculations within a reasonable 
scope. (For instance, we have selected a system size un- 
der 100,000 orbital sites to facilitate replication of our re- 
sults). Since trigonal warping becomes significant above 
E = 0.4t, we have selected the energy of 0.2£ for all states 
in our analysis to maximize the number of wavelengths 
within a small graphene system while maintaining the 




Figure 7: Low-energy graphene states require additional tools 
to fully grasp the classical dynamics. The Husimi map for the 
if'-valley is plotted for four eigenstates of a closed circular 
system with 71934 orbital sites at energies around E = 0.2t. 
All three calculations use coherent states with relative un- 
certainty Ak/k = 20%, with breadth indicated by the dou- 
ble arrows on the right. From left-to-right: the Husimi flux, 
multi-modal analysis, and the wavefunction. The divergence 
of the Husimi flux is indicated in green (red) for positive (neg- 
ative) values. In blue, the angular deflection. 



same physics from energies closer to the Dirac point. 

Fig. [7] shows four eigenstates of the circular graphene 
flake. Like the free-particle circular well, eigenstates of 
the graphene circular flake resemble eigenstates of the 
angular momentum operator (see Mason et al. [24J for di- 
rect comparisons and Husimi maps). For instance, the 
wavefunctions in Figs. [7k- b are radial-dominant, while 
the wavefunction in Fig7[7]i is angular-dominant. These 
observations carry over to the dynamics of the wavefunc- 
tions revealed by the multi-modal analysis for the K' 
valley, which shows radially-oriented paths in Figs. 7i- 
b and circular paths skimming the boundary in Fig. 71. 
Fig. [T]: shows a state with a mixture of radial and angular 
components; in the multi-modal analysis, this appears as 
straight paths between boundary points highlighted by 
the angular deflection. 

Unlike free-particle circular wells, however, the lat- 
tice sampling on the honeycomb lattice breaks circular 
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Figure 8: In parts (a) and (b), the same information is plotted 
as in Fig. [7] but for the Wimmer system (see Fig. [6}, with 
96425 orbital sites. These states also have energies near E — 
0.2£ and are represented by coherent states of uncertainty 
Ak/k — 20% with breadth indicated by the double arrows. 



symmetry and replaces it with six-fold symmetry. Be- 
cause eigenstates of the system emphasize certain bound- 
ary conditions, the manner in which each state estab- 
lishes itself strongly varies. For instance, the two radial- 
dominant states in Figs.[7^i-b exhibit intravalley (a) or in- 
tervalley (b) scattering. Accordingly, the locations where 
the rays terminate on the boundary correlate with zig-zag 
and armchair boundaries respectively. The wider spread 
in angular deflection in Fig. [T^i corroborates Akhmerov 
and Beenakker[34j, showing that that intravalley scatter- 
ing occurs over a larger set of boundaries than intervalley 
scattering. 

Because each valley reflects back to itself in Fig. [T^i, 
there is no net flow of either valley in the bulk of the 
system. As a result, the multi- modal analysis shows 
counter-propagating flows, and the Husimi flux (Eq.[5| is 
zero except at the center, where slight offsets in trajecto- 
ries form characteristic vortices. In Fig. [7]3, on the other 
hand, each ray in the wavefunction is associated with a 
distinct source and drain, which is evident in both the 
multi-modal analysis and the Husimi flux. 

In Figs. [7]>d, the locations of sources and drains for 
the K' valley are reversed from Fig. [7)3. However, the 
roles that inter-valley scattering play in these states is 
less clear; rather, inter- and intra- valley scattering domi- 
nate these wavefunctions. In Fig. [7]:, this can been by the 
emphasis of angular deflection along the zig-zag bound- 
aries, which do not show any divergence. In Fig.[7]i, even 
though the wavefunction and the multi-modal analysis 
clearly emphasize a classical path that skims the bound- 
ary, the path actually flips between each valley each time 
it encounters an inter- valley scatterer. For both states, 
the various trajectories merge to form vortices in the 
Husimi flux, with sources and drains at armchair edges. 

When the circular flake is distorted, as in the Wimmer 
system (Figs. [6] and [8]), inter- and intra- valley scatterers 
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Figure 9: An extremely small "rooftop" graphene flake at en- 
ergy E — 0.0015735£ showing two edge states at the top and 
bottom boundaries which tunnel into each other. At top, the 
full wavefunction, at middle, divergence is indicated in green 
and red, and a schematic of the Husimi flux for the K' val- 
ley is shown. At bottom, the Fourier transform of the state 
is shown, with the contour line used to generate the Husimi 
map in white, set at an arbitrary energy in order to maximize 
the intersection of the contour with the Fourier-transform am- 
plitude. The spread of the wavepacket used to generate this 
map is indicated by the double arrows. 



are re-arranged and re-sized as a function of the local 
radius of curvature of the boundary. 

Figs.[8^-b show two eigenstates of the Wimmer system. 
The boundary conditions for these states most closely 
resemble Fig. |5j since sources and drains appear next to 
each other. This is a signature of mixed scattering - both 
inter- and intra- valley scattering occur in various propor- 
tions at these points. For example, the multi- modal anal- 
ysis in Fig. [8^1 shows a triangular path, but not all legs of 
the triangle are equally strong, corresponding to various 
degrees of absorption and reflection at each scattering 
point which can be seen in the divergence. 

Edge-states are a set of zero-energy surface states 
that are strongly localized to zig-zag boundaries and 
potentially long-lived [9j, and since they can be used 
as modes of transport [TUl [12] and be strongly spin- 
polarized|TTJ [13], they have been proposed a candidate 
for spin-tronics devices [9HT3j. But because edge states 
exhibit a different dispersion relation than the two val- 
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Figure 10: System properties of the scattering density matrix 
p around the Fano resonance centered at E — 1.9582eV for 
the open system in the inset. Top: The transmission profile 
across the two leads, with the closed- system eigenstate en- 
ergy at E = 1.9579eV, corresponding to the eigenstate at in- 
dex 1483 (below), indicated by the vertical grey line. Middle: 
Diagonalizing the density matrix produces a handful of non- 
trivial scattering wavefunctions in its eigenvectors. The eigen- 
values of these vectors, which correspond to their measure- 
ment probability, are graphed. The wavefunction associated 
with the closed-system eigenstate hybridizing with the direct 
channel peaks strongly around the Fano resonance. Bottom: 
The density matrix is projected onto the closed-system eigen- 
states, showing that eigenstate 1483 strongly peaks at the 
Fano resonance. 



leys in the bulk, they cannot be "sensed" by the K' or 
K valley Husimi projections. Instead, the Husimi map 
can be generated using wavevectors appropriate to the 
edge states, which shows them as standing waves on the 
surface (see Fig. [9|. As noted in Wimmer et al. [35], it 
is possible for edge states to tunnel into each other us- 
ing bulk states as a medium, but we have found that 
K' or K valley Husimi maps of bulk states which hy- 
bridize with them are indistinguishable from their non- 
hybridized counterparts 



C. Fano Resonance 

This section addresses Fano resonance [28J in graphene 
systems, a conductance phenomenon that occurs as a 
result of interference between a direct state (conduc- 
tance channel) and a quasi-bound indirect state similar 
to the eigenstates this paper has examined. Fano reso- 
nances are an ideal case study for the Husimi map, not 



only because they are ubiquitous in theory [36] [37] and 
experiments [38 -40J, but also because their behavior is 
well-understood [30, 41-46J. However, Fano resonances 
in graphene quantum dots are less well characterized [47- 
150] and lack a comprehensive theory relating boundary 
conditions to bulk state behavior in graphene. 

To study Fano resonance, we first compute a scattering 
wavefunction using the recursive numerical Green's func- 
tion method described in Mason et al. [51]. This method 
produces a scattering density matrix p, which is diagonal- 
ized. Each eigenvector corresponds to a scattering wave- 
function, which has an associated eigenvalue indicating 
its measurement probability (Fig. [ToJ middle). We focus 
on the resonance study in Huang et al.[27\. 

The resonance in Fig. [10] is associated with the eigen- 
state from Fig.[5]of the closed billiard system. This eigen- 
state couples only weakly to leads which are attached at 
its sides (shown in the inset of Fig. 10). This makes 



it possible for a scattering electron to enter the system 
through a direct channel but then become trapped in a 
quasi-bound state related to the eigenstate, causing the 
density of states projected onto the eigenstate to strongly 
peak near its eigenenergy (Fig. [ToJ bottom) . As the sys- 
tem energy sweeps across the eigenenergy, the phase of 
the eigenstate component shifts through 7r, causing it to 
interfere negatively and then positively with the direct 
channel, giving rise to the distinctive Fano curve (Fig. 10 , 
top). As a result, the scattering wavefunction with the 
largest measurement probability is in fact a hybridized 
state between the closed-system eigenstate and the di- 
rect channel, and its probability peaks around an energy 
near, but not exactly the same as, the eigenstate energy 
(Fig. [ToJ middle) . The shift in energy arises as a pertur- 
bation from the leads. 

For closed graphene systems, the two valleys satisfy 
time-reversal symmetry as an analytical consequence of 
lattice sampling on the honeycomb lattice. As a result, 
trajectories in one valley are exactly reversed from the 
other valley, in analogy with free-particle systems where 
opposing trajectories cancel each other produce zero flux. 
This observation allows us to remove the time-reversal 
symmetry of a scattering wavefunction by summing the 
projections for both valleys, revealing the time-reversal 
asymmetric part of the wavefunction. 



Fig. 11 shows the results of adding the Husimi flux 
maps of both valleys at two energies, below and above 
resonance. We find sources and drains in the summed 
Husimi flux map at the corners of the system where the 
classical paths of the if '-valley Husimi map (Fig. [5]) re- 
flect off the system boundary. 

To understand why, we consider that during trans- 
mission, quasiparticles enter from the left incoming lead 
and exit through the right outgoing lead. However, near 
resonance, the wavefunction is strongly weighted by the 
closed-system eigenstate, which has no net quasiparticle 
current. Husimi maps for either valley also reflect this 
fact: they are indistinguishable from the Husimi maps of 
the closed-system eigenstate in Fig.[5| and the two valleys 



8 




1 I 

v \ \ \ 



miw 



Figure 11: Above and below the Fano resonance in Figs. [l0| 
(inset), the time-reversal symmetry between the K and K 
valleys is lifted, making it possible to add the Husimi flux for 
both valleys to measure valley-polarized current. Above, the 
Husimi flux maps of both valleys are added for the scatter- 
ing wavefunction at energies E — 1.9582£ and 1.9586£, with 
Ak/k = 30%. Below, the probability flux, convolved with a 
Gaussian kernel of the same size as the coherent state. At 
energies this close to resonance, the wavefunction does not 
visually change from the closed-system eigenstate in the in- 
set, but the residual current that occurs near these resonances 
switches direction across resonance. 



are inverse images of each other. 

But the Husimi maps for the two valleys don't ex- 
actly cancel each other out. When we add them together 
to reveal the time-reversal asymmetric behavior of the 
wavefunction, the residual shows sources and drains of 
net quasiparticle flow which are strongly related to the 
Husimi maps for each valley, and do not show left-to- 
right transmission. Instead, the summed Husimi flux 
map shows the influence of transmission on the strongly- 
emphasized classical paths underlying the closed-system 
eigenstate. 

To compare the summed Husimi flux map to the tra- 
ditional flux, we consider the probability flow between 
two adjacent carbon atom sites called the bond current, 
defined as 



where Hij and Gfj(E) are the off-diagonal components 
of the Hamiltonian and the electron correlation function 
between orbital sites i and j \%T\ [5^] . The electron cor- 
relation function is proportional to the density matrix, 
but in our calculations, we examine just one scattering 
state, so that Gfj oc where ipi is the scattering state 
probability amplitude at orbital site i. We can obtain a 
finite-difference analog of the continuum flux operator by 
defining 



-^3 | 



|2 ' 



(16) 



which computes the vector sum of each bond current as- 
sociated with a given orbital[53j. 



4e 



Im [HyCfyiE)], 



(15) 



Convolving the flux defined in Eq. 16 with a Gaussian 
kernel of the same spread as the coherent state used to 
generate the Husimi map creates an analog to the Husimi 
flux, except that the convolved flux does not distinguish 
among valleys. We show the convolved flux at the bottom 
of Fig. [TT] and find that it forms vortices which correlates 
with the summed Husimi flux maps, and also fails to show 
the left-to-right flow responsible for transmission. 

This behavior is directly analogous to flux in contin- 
uum systems, where flux vortices above and below res- 
onance show local variations of flow but not the left-to- 
right drift velocity responsible for transmission. We can 
recover the left-to-right flow only by examining the sys- 
tem at larger scales using larger Gaussian spreads (not 
shown) [24J. Because of the tt phase shift of the indi- 
rect channel across resonance, local flows reverse direc- 
tion above and below resonance, but they do not affect 
the left-to-right flow at larger scales except exactly on 
resonance. 

The stable orbits that underly the indirect channel, 
shown in Fig. [5j can be dramatically disturbed by slight 
modifications of the boundary where the classical paths 
reflect off the boundary. The original authors Huang et 
al.[27\ examined the relationship between system symme- 
try and strength of the Fano resonances by slightly mod- 
ifying the system boundary at the black circle in Fig. [5j 
and demonstrated that some resonances were drastically 
reduced by this modification. 

We have chosen the resonance in this study because 
the Fano resonance profile associated with it was among 
the most-reduced as a result of their system modifica- 
tion, and our analysis provides a clear picture as to why: 
the system is perturbed precisely at the boundary where 
the eigenstate in Fig. [5] has the largest probability am- 
plitude. With the semiclassical picture, we are able to 
add to this finding an intuitive understanding: by dis- 
turbing the reflection angle at the exact point where the 
two valleys scatter, each time the electron scatters off 
that point some of its probability leaves the stable orbit. 
The authors effectively introduced a leak into the orbit, 
reducing its lifetime and the strength of its resonance 
considerably. 
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IV. CONCLUSIONS 

We have examined the semiclassical behavior of 
graphene systems using a generalized technique that 
produces a vector field from projections onto coherent 
states, forming an infinitely tunable bridge between the 
large-scale Dirac effective field theory and the underlying 
atomistic model |9J. We have used this technique, called 
the Husimi map, to examine the relationship between 
graphene boundary types and the classical dynamics of 
quasiparticles in each valley of the honeycomb dispersion 
relation, looking at states with energies both close to and 
far from the Dirac point. We have shown that closed- 
system eigenstates are associated with valley-polarized 
currents with zero net quasiparticle production. We show 
that Fano resonance are associated with an asymmetri- 



cal flow of quasiparticles strongly related to the valley- 
polarized currents of closed- system states, which has im- 
plications for applications in "valley tronic" devices [54J. 
The ubiquity of this phenomenon in the systems we have 
studied suggests that they could appear in future exper- 
iments, and provides a motivation for further theoretical 
and experimental work. 
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